home *** CD-ROM | disk | FTP | other *** search
/ Games of Daze / Infomagic - Games of Daze (Summer 1995) (Disc 1 of 2).iso / x2ftp / msdos / source / gemsiv / implicit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-05-16  |  24.0 KB  |  816 lines
  1. /*
  2.  * C code from the article
  3.  * "An Implicit Surface Polygonizer"
  4.  * by Jules Bloomenthal, jbloom@beauty.gmu.edu
  5.  * in "Graphics Gems IV", Academic Press, 1994
  6.  */
  7.  
  8. /* implicit.c
  9.  *     an implicit surface polygonizer, translated from Mesa
  10.  *     applications should call polygonize()
  11.  *
  12.  * to compile a test program for ASCII output:
  13.  *     cc implicit.c -o implicit -lm
  14.  *
  15.  * to compile a test program for display on an SGI workstation:
  16.  *     cc -DSGIGFX implicit.c -o implicit -lgl_s -lm
  17.  *
  18.  * Authored by Jules Bloomenthal, Xerox PARC.
  19.  * Copyright (c) Xerox Corporation, 1991.  All rights reserved.
  20.  * Permission is granted to reproduce, use and distribute this code for
  21.  * any and all purposes, provided that this notice appears in all copies. */
  22.  
  23. #include <stdlib.h>
  24. #include <math.h>
  25. #include <stdio.h>
  26. #include <sys/types.h>
  27.  
  28. #define TET    0  /* use tetrahedral decomposition */
  29. #define NOTET    1  /* no tetrahedral decomposition  */
  30.  
  31. #define RES    10 /* # converge iterations    */
  32.  
  33. #define L    0  /* left direction:    -x, -i */
  34. #define R    1  /* right direction:    +x, +i */
  35. #define B    2  /* bottom direction: -y, -j */
  36. #define T    3  /* top direction:    +y, +j */
  37. #define N    4  /* near direction:    -z, -k */
  38. #define F    5  /* far direction:    +z, +k */
  39. #define LBN    0  /* left bottom near corner  */
  40. #define LBF    1  /* left bottom far corner   */
  41. #define LTN    2  /* left top near corner     */
  42. #define LTF    3  /* left top far corner      */
  43. #define RBN    4  /* right bottom near corner */
  44. #define RBF    5  /* right bottom far corner  */
  45. #define RTN    6  /* right top near corner    */
  46. #define RTF    7  /* right top far corner     */
  47.  
  48. /* the LBN corner of cube (i, j, k), corresponds with location
  49.  * (start.x+(i-.5)*size, start.y+(j-.5)*size, start.z+(k-.5)*size) */
  50.  
  51. #define RAND()        ((rand()&32767)/32767.)    /* random number between 0 and 1 */
  52. #define HASHBIT        (5)
  53. #define HASHSIZE    (size_t)(1<<(3*HASHBIT))   /* hash table size (32768) */
  54. #define MASK        ((1<<HASHBIT)-1)
  55. #define HASH(i,j,k) ((((((i)&MASK)<<HASHBIT)|((j)&MASK))<<HASHBIT)|((k)&MASK))
  56. #define BIT(i, bit) (((i)>>(bit))&1)
  57. #define FLIP(i,bit) ((i)^1<<(bit)) /* flip the given bit of i */
  58.  
  59. typedef struct point {           /* a three-dimensional point */
  60.     double x, y, z;           /* its coordinates */
  61. } POINT;
  62.  
  63. typedef struct test {           /* test the function for a signed value */
  64.     POINT p;               /* location of test */
  65.     double value;           /* function value at p */
  66.     int ok;               /* if value is of correct sign */
  67. } TEST;
  68.  
  69. typedef struct vertex {           /* surface vertex */
  70.     POINT position, normal;       /* position and surface normal */
  71. } VERTEX;
  72.  
  73. typedef struct vertices {       /* list of vertices in polygonization */
  74.     int count, max;           /* # vertices, max # allowed */
  75.     VERTEX *ptr;           /* dynamically allocated */
  76. } VERTICES;
  77.  
  78. typedef struct corner {           /* corner of a cube */
  79.     int i, j, k;           /* (i, j, k) is index within lattice */
  80.     double x, y, z, value;       /* location and function value */
  81. } CORNER;
  82.  
  83. typedef struct cube {           /* partitioning cell (cube) */
  84.     int i, j, k;           /* lattice location of cube */
  85.     CORNER *corners[8];           /* eight corners */
  86. } CUBE;
  87.  
  88. typedef struct cubes {           /* linked list of cubes acting as stack */
  89.     CUBE cube;               /* a single cube */
  90.     struct cubes *next;           /* remaining elements */
  91. } CUBES;
  92.  
  93. typedef struct centerlist {       /* list of cube locations */
  94.     int i, j, k;           /* cube location */
  95.     struct centerlist *next;       /* remaining elements */
  96. } CENTERLIST;
  97.  
  98. typedef struct cornerlist {       /* list of corners */
  99.     int i, j, k;           /* corner id */
  100.     double value;           /* corner value */
  101.     struct cornerlist *next;       /* remaining elements */
  102. } CORNERLIST;
  103.  
  104. typedef struct edgelist {       /* list of edges */
  105.     int i1, j1, k1, i2, j2, k2;       /* edge corner ids */
  106.     int vid;               /* vertex id */
  107.     struct edgelist *next;       /* remaining elements */
  108. } EDGELIST;
  109.  
  110. typedef struct intlist {       /* list of integers */
  111.     int i;               /* an integer */
  112.     struct intlist *next;       /* remaining elements */
  113. } INTLIST;
  114.  
  115. typedef struct intlists {       /* list of list of integers */
  116.     INTLIST *list;           /* a list of integers */
  117.     struct intlists *next;       /* remaining elements */
  118. } INTLISTS;
  119.  
  120. typedef struct process {       /* parameters, function, storage */
  121.     double (*function)();       /* implicit surface function */
  122.     int (*triproc)();           /* triangle output function */
  123.     double size, delta;           /* cube size, normal delta */
  124.     int bounds;               /* cube range within lattice */
  125.     POINT start;           /* start point on surface */
  126.     CUBES *cubes;           /* active cubes */
  127.     VERTICES vertices;           /* surface vertices */
  128.     CENTERLIST **centers;       /* cube center hash table */
  129.     CORNERLIST **corners;       /* corner value hash table */
  130.     EDGELIST **edges;           /* edge and vertex id hash table */
  131. } PROCESS;
  132.  
  133. void *calloc();
  134. char *mycalloc();
  135.  
  136.  
  137. /**** A Test Program ****/
  138.  
  139.  
  140. /* torus: a torus with major, minor radii = 0.5, 0.1, try size = .05 */
  141.  
  142. double torus (x, y, z)
  143. double x, y, z;
  144. {
  145.     double x2 = x*x, y2 = y*y, z2 = z*z;
  146.     double a = x2+y2+z2+(0.5*0.5)-(0.1*0.1);
  147.     return a*a-4.0*(0.5*0.5)*(y2+z2);
  148. }
  149.  
  150.  
  151. /* sphere: an inverse square function (always positive) */
  152.  
  153. double sphere (x, y, z)
  154. double x, y, z;
  155. {
  156.     double rsq = x*x+y*y+z*z;
  157.     return 1.0/(rsq < 0.00001? 0.00001 : rsq);
  158. }
  159.  
  160.  
  161. /* blob: a three-pole blend function, try size = .1 */
  162.  
  163. double blob (x, y, z)
  164. double x, y, z;
  165. {
  166.     return 4.0-sphere(x+1.0,y,z)-sphere(x,y+1.0,z)-sphere(x,y,z+1.0);
  167. }
  168.  
  169.  
  170. #ifdef SGIGFX /**************************************************************/
  171.  
  172. #include "gl.h"
  173.  
  174. /* triangle: called by polygonize() for each triangle; set SGI lines */
  175.  
  176. triangle (i1, i2, i3, vertices)
  177. int i1, i2, i3;
  178. VERTICES vertices;
  179. {
  180.     float v[3];
  181.     int i, ids[3];
  182.     ids[0] = i1;
  183.     ids[1] = i2;
  184.     ids[2] = i3;
  185.     bgnclosedline();
  186.     for (i = 0; i < 3; i++) {
  187.     POINT *p = &vertices.ptr[ids[i]].position;
  188.     v[0] = p->x; v[1] = p->y; v[2] = p->z;
  189.     v3f(v);
  190.     }
  191.     endclosedline();
  192.     return 1;
  193. }
  194.  
  195.  
  196. /* main: call polygonize() with torus function
  197.  * display lines on SGI */
  198.  
  199. main ()
  200. {
  201.     char *err, *polygonize();
  202.  
  203.     keepaspect(1, 1);
  204.     winopen("implicit");
  205.     doublebuffer();
  206.     gconfig();
  207.     perspective(450, 1.0/1.0, 0.1, 10.0);
  208.     color(7);
  209.     clear();
  210.     swapbuffers();
  211.     makeobj(1);
  212.     if ((err = polygonize(torus, .05, 20, 0.,0.,0., triangle, TET)) != NULL) {
  213.     fprintf(stderr, "%s\n", err);
  214.     exit(1);
  215.     }
  216.     closeobj();
  217.     translate(0.0, 0.0, -2.0);
  218.     pushmatrix();
  219.     while(1) { /* spin the object */
  220.     reshapeviewport();
  221.     color(7);
  222.     clear();
  223.     color(0);
  224.     callobj(1);
  225.     rot(0.8, 'x');
  226.     rot(0.3, 'y');
  227.     rot(0.1, 'z');
  228.     swapbuffers();
  229.  
  230.     }
  231. }
  232.  
  233. #else /***********************************************************************/
  234.  
  235. int gntris;         /* global needed by application */
  236. VERTICES gvertices;  /* global needed by application */
  237.  
  238.  
  239. /* triangle: called by polygonize() for each triangle; write to stdout */
  240.  
  241. triangle (i1, i2, i3, vertices)
  242. int i1, i2, i3;
  243. VERTICES vertices;
  244. {
  245.     gvertices = vertices;
  246.     gntris++;
  247.     fprintf(stdout, "%d %d %d\n", i1, i2, i3);
  248.     return 1;
  249. }
  250.  
  251.  
  252. /* main: call polygonize() with torus function
  253.  * write points-polygon formatted data to stdout */
  254.  
  255. main ()
  256.     {
  257.     int i;
  258.     char *err, *polygonize();
  259.     gntris = 0;
  260.     fprintf(stdout, "triangles\n\n");
  261.     if ((err = polygonize(torus, .05, 20, 0.,0.,0., triangle, TET)) != NULL) {
  262.     fprintf(stdout, "%s\n", err);
  263.     exit(1);
  264.     }
  265.     fprintf(stdout, "\n%d triangles, %d vertices\n", gntris, gvertices.count);
  266.     fprintf(stdout, "\nvertices\n\n");
  267.     for (i = 0; i < gvertices.count; i++) {
  268.     VERTEX v;
  269.     v = gvertices.ptr[i];
  270.     fprintf(stdout, "%f  %f     %f\t%f     %f  %f\n",
  271.         v.position.x, v.position.y,     v.position.z,
  272.         v.normal.x,      v.normal.y,     v.normal.z);
  273.     }
  274.     fprintf(stderr, "%d triangles, %d vertices\n", gntris, gvertices.count);
  275.     exit(0);
  276. }
  277.  
  278. #endif /**********************************************************************/
  279.  
  280.  
  281. /**** An Implicit Surface Polygonizer ****/
  282.  
  283.  
  284. /* polygonize: polygonize the implicit surface function
  285.  *   arguments are:
  286.  *     double function (x, y, z)
  287.  *         double x, y, z (an arbitrary 3D point)
  288.  *         the implicit surface function
  289.  *         return negative for inside, positive for outside
  290.  *     double size
  291.  *         width of the partitioning cube
  292.  *     int bounds
  293.  *         max. range of cubes (+/- on the three axes) from first cube
  294.  *     double x, y, z
  295.  *         coordinates of a starting point on or near the surface
  296.  *         may be defaulted to 0., 0., 0.
  297.  *     int triproc (i1, i2, i3, vertices)
  298.  *         int i1, i2, i3 (indices into the vertex array)
  299.  *         VERTICES vertices (the vertex array, indexed from 0)
  300.  *         called for each triangle
  301.  *         the triangle coordinates are (for i = i1, i2, i3):
  302.  *         vertices.ptr[i].position.x, .y, and .z
  303.  *         vertices are ccw when viewed from the out (positive) side
  304.  *         in a left-handed coordinate system
  305.  *         vertex normals point outwards
  306.  *         return 1 to continue, 0 to abort
  307.  *     int mode
  308.  *         TET: decompose cube and polygonize six tetrahedra
  309.  *         NOTET: polygonize cube directly
  310.  *   returns error or NULL
  311.  */
  312.  
  313. char *polygonize (function, size, bounds, x, y, z, triproc, mode)
  314. double (*function)(), size, x, y, z;
  315. int bounds, (*triproc)(), mode;
  316. {
  317.     PROCESS p;
  318.     int n, noabort;
  319.     CORNER *setcorner();
  320.     TEST in, out, find();
  321.  
  322.     p.function = function;
  323.     p.triproc = triproc;
  324.     p.size = size;
  325.     p.bounds = bounds;
  326.     p.delta = size/(double)(RES*RES);
  327.  
  328.     /* allocate hash tables and build cube polygon table: */
  329.     p.centers = (CENTERLIST **) mycalloc(HASHSIZE,sizeof(CENTERLIST *));
  330.     p.corners = (CORNERLIST **) mycalloc(HASHSIZE,sizeof(CORNERLIST *));
  331.     p.edges =    (EDGELIST   **) mycalloc(2*HASHSIZE,sizeof(EDGELIST *));
  332.     makecubetable();
  333.  
  334.     /* find point on surface, beginning search at (x, y, z): */
  335.     srand(1);
  336.     in = find(1, &p, x, y, z);
  337.     out = find(0, &p, x, y, z);
  338.     if (!in.ok || !out.ok) return "can't find starting point";
  339.     converge(&in.p, &out.p, in.value, p.function, &p.start);
  340.  
  341.     /* push initial cube on stack: */
  342.     p.cubes = (CUBES *) mycalloc(1, sizeof(CUBES)); /* list of 1 */
  343.     p.cubes->cube.i = p.cubes->cube.j = p.cubes->cube.k = 0;
  344.     p.cubes->next = NULL;
  345.  
  346.     /* set corners of initial cube: */
  347.     for (n = 0; n < 8; n++)
  348.     p.cubes->cube.corners[n] = setcorner(&p, BIT(n,2), BIT(n,1), BIT(n,0));
  349.  
  350.     p.vertices.count = p.vertices.max = 0; /* no vertices yet */
  351.     p.vertices.ptr = NULL;
  352.  
  353.     setcenter(p.centers, 0, 0, 0);
  354.  
  355.     while (p.cubes != NULL) { /* process active cubes till none left */
  356.     CUBE c;
  357.     CUBES *temp = p.cubes;
  358.     c = p.cubes->cube;
  359.  
  360.     noabort = mode == TET?
  361.            /* either decompose into tetrahedra and polygonize: */
  362.            dotet(&c, LBN, LTN, RBN, LBF, &p) &&
  363.            dotet(&c, RTN, LTN, LBF, RBN, &p) &&
  364.            dotet(&c, RTN, LTN, LTF, LBF, &p) &&
  365.            dotet(&c, RTN, RBN, LBF, RBF, &p) &&
  366.            dotet(&c, RTN, LBF, LTF, RBF, &p) &&
  367.            dotet(&c, RTN, LTF, RTF, RBF, &p)
  368.            :
  369.            /* or polygonize the cube directly: */
  370.            docube(&c, &p);
  371.     if (! noabort) return "aborted";
  372.  
  373.     /* pop current cube from stack */
  374.     p.cubes = p.cubes->next;
  375.     free((char *) temp);
  376.     /* test six face directions, maybe add to stack: */
  377.     testface(c.i-1, c.j, c.k, &c, L, LBN, LBF, LTN, LTF, &p);
  378.     testface(c.i+1, c.j, c.k, &c, R, RBN, RBF, RTN, RTF, &p);
  379.     testface(c.i, c.j-1, c.k, &c, B, LBN, LBF, RBN, RBF, &p);
  380.     testface(c.i, c.j+1, c.k, &c, T, LTN, LTF, RTN, RTF, &p);
  381.     testface(c.i, c.j, c.k-1, &c, N, LBN, LTN, RBN, RTN, &p);
  382.     testface(c.i, c.j, c.k+1, &c, F, LBF, LTF, RBF, RTF, &p);
  383.     }
  384.     return NULL;
  385. }
  386.  
  387.  
  388. /* testface: given cube at lattice (i, j, k), and four corners of face,
  389.  * if surface crosses face, compute other four corners of adjacent cube
  390.  * and add new cube to cube stack */
  391.  
  392. testface (i, j, k, old, face, c1, c2, c3, c4, p)
  393. CUBE *old;
  394. PROCESS *p;
  395. int i, j, k, face, c1, c2, c3, c4;
  396. {
  397.     CUBE new;
  398.     CUBES *oldcubes = p->cubes;
  399.     CORNER *setcorner();
  400.     static int facebit[6] = {2, 2, 1, 1, 0, 0};
  401.     int n, pos = old->corners[c1]->value > 0.0 ? 1 : 0, bit = facebit[face];
  402.  
  403.     /* test if no surface crossing, cube out of bounds, or already visited: */
  404.     if ((old->corners[c2]->value > 0) == pos &&
  405.     (old->corners[c3]->value > 0) == pos &&
  406.     (old->corners[c4]->value > 0) == pos) return;
  407.     if (abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;
  408.     if (setcenter(p->centers, i, j, k)) return;
  409.  
  410.     /* create new cube: */
  411.     new.i = i;
  412.     new.j = j;
  413.     new.k = k;
  414.     for (n = 0; n < 8; n++) new.corners[n] = NULL;
  415.     new.corners[FLIP(c1, bit)] = old->corners[c1];
  416.     new.corners[FLIP(c2, bit)] = old->corners[c2];
  417.     new.corners[FLIP(c3, bit)] = old->corners[c3];
  418.     new.corners[FLIP(c4, bit)] = old->corners[c4];
  419.     for (n = 0; n < 8; n++)
  420.     if (new.corners[n] == NULL)
  421.         new.corners[n] = setcorner(p, i+BIT(n,2), j+BIT(n,1), k+BIT(n,0));
  422.  
  423.     /*add cube to top of stack: */
  424.     p->cubes = (CUBES *) mycalloc(1, sizeof(CUBES));
  425.     p->cubes->cube = new;
  426.     p->cubes->next = oldcubes;
  427. }
  428.  
  429.  
  430. /* setcorner: return corner with the given lattice location
  431.    set (and cache) its function value */
  432.  
  433. CORNER *setcorner (p, i, j, k)
  434. int i, j, k;
  435. PROCESS *p;
  436. {
  437.     /* for speed, do corner value caching here */
  438.     CORNER *c = (CORNER *) mycalloc(1, sizeof(CORNER));
  439.     int index = HASH(i, j, k);
  440.     CORNERLIST *l = p->corners[index];
  441.     c->i = i; c->x = p->start.x+((double)i-.5)*p->size;
  442.     c->j = j; c->y = p->start.y+((double)j-.5)*p->size;
  443.     c->k = k; c->z = p->start.z+((double)k-.5)*p->size;
  444.     for (; l != NULL; l = l->next)
  445.     if (l->i == i && l->j == j && l->k == k) {
  446.         c->value = l->value;
  447.         return c;
  448.         }
  449.     l = (CORNERLIST *) mycalloc(1, sizeof(CORNERLIST));
  450.     l->i = i; l->j = j; l->k = k;
  451.     l->value = c->value = p->function(c->x, c->y, c->z);
  452.     l->next = p->corners[index];
  453.     p->corners[index] = l;
  454.     return c;
  455. }
  456.  
  457.  
  458. /* find: search for point with value of given sign (0: neg, 1: pos) */
  459.  
  460. TEST find (sign, p, x, y, z)
  461. int sign;
  462. PROCESS *p;
  463. double x, y, z;
  464. {
  465.     int i;
  466.     TEST test;
  467.     double range = p->size;
  468.     test.ok = 1;
  469.     for (i = 0; i < 10000; i++) {
  470.     test.p.x = x+range*(RAND()-0.5);
  471.     test.p.y = y+range*(RAND()-0.5);
  472.     test.p.z = z+range*(RAND()-0.5);
  473.     test.value = p->function(test.p.x, test.p.y, test.p.z);
  474.     if (sign == (test.value > 0.0)) return test;
  475.     range = range*1.0005; /* slowly expand search outwards */
  476.     }
  477.     test.ok = 0;
  478.     return test;
  479. }
  480.  
  481.  
  482. /**** Tetrahedral Polygonization ****/
  483.  
  484.  
  485. /* dotet: triangulate the tetrahedron
  486.  * b, c, d should appear clockwise when viewed from a
  487.  * return 0 if client aborts, 1 otherwise */
  488.  
  489. int dotet (cube, c1, c2, c3, c4, p)
  490. CUBE *cube;
  491. int c1, c2, c3, c4;
  492. PROCESS *p;
  493. {
  494.     CORNER *a = cube->corners[c1];
  495.     CORNER *b = cube->corners[c2];
  496.     CORNER *c = cube->corners[c3];
  497.     CORNER *d = cube->corners[c4];
  498.     int index = 0, apos, bpos, cpos, dpos, e1, e2, e3, e4, e5, e6;
  499.     if (apos = (a->value > 0.0)) index += 8;
  500.     if (bpos = (b->value > 0.0)) index += 4;
  501.     if (cpos = (c->value > 0.0)) index += 2;
  502.     if (dpos = (d->value > 0.0)) index += 1;
  503.     /* index is now 4-bit number representing one of the 16 possible cases */
  504.     if (apos != bpos) e1 = vertid(a, b, p);
  505.     if (apos != cpos) e2 = vertid(a, c, p);
  506.     if (apos != dpos) e3 = vertid(a, d, p);
  507.     if (bpos != cpos) e4 = vertid(b, c, p);
  508.     if (bpos != dpos) e5 = vertid(b, d, p);
  509.     if (cpos != dpos) e6 = vertid(c, d, p);
  510.     /* 14 productive tetrahedral cases (0000 and 1111 do not yield polygons */
  511.     switch (index) {
  512.     case 1:     return p->triproc(e5, e6, e3, p->vertices);
  513.     case 2:     return p->triproc(e2, e6, e4, p->vertices);
  514.     case 3:     return p->triproc(e3, e5, e4, p->vertices) &&
  515.             p->triproc(e3, e4, e2, p->vertices);
  516.     case 4:     return p->triproc(e1, e4, e5, p->vertices);
  517.     case 5:     return p->triproc(e3, e1, e4, p->vertices) &&
  518.             p->triproc(e3, e4, e6, p->vertices);
  519.     case 6:     return p->triproc(e1, e2, e6, p->vertices) &&
  520.             p->triproc(e1, e6, e5, p->vertices);
  521.     case 7:     return p->triproc(e1, e2, e3, p->vertices);
  522.     case 8:     return p->triproc(e1, e3, e2, p->vertices);
  523.     case 9:     return p->triproc(e1, e5, e6, p->vertices) &&
  524.             p->triproc(e1, e6, e2, p->vertices);
  525.     case 10: return p->triproc(e1, e3, e6, p->vertices) &&
  526.             p->triproc(e1, e6, e4, p->vertices);
  527.     case 11: return p->triproc(e1, e5, e4, p->vertices);
  528.     case 12: return p->triproc(e3, e2, e4, p->vertices) &&
  529.             p->triproc(e3, e4, e5, p->vertices);
  530.     case 13: return p->triproc(e6, e2, e4, p->vertices);
  531.     case 14: return p->triproc(e5, e3, e6, p->vertices);
  532.     }
  533.     return 1;
  534. }
  535.  
  536.  
  537. /**** Cubical Polygonization (optional) ****/
  538.  
  539.  
  540. #define LB    0  /* left bottom edge    */
  541. #define LT    1  /* left top edge    */
  542. #define LN    2  /* left near edge    */
  543. #define LF    3  /* left far edge    */
  544. #define RB    4  /* right bottom edge */
  545. #define RT    5  /* right top edge    */
  546. #define RN    6  /* right near edge    */
  547. #define RF    7  /* right far edge    */
  548. #define BN    8  /* bottom near edge    */
  549. #define BF    9  /* bottom far edge    */
  550. #define TN    10 /* top near edge    */
  551. #define TF    11 /* top far edge    */
  552.  
  553. static INTLISTS *cubetable[256];
  554.  
  555. /*            edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
  556. static int corner1[12]       = {LBN,LTN,LBN,LBF,RBN,RTN,RBN,RBF,LBN,LBF,LTN,LTF};
  557. static int corner2[12]       = {LBF,LTF,LTN,LTF,RBF,RTF,RTN,RTF,RBN,RBF,RTN,RTF};
  558. static int leftface[12]       = {B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F};
  559.                  /* face on left when going corner1 to corner2 */
  560. static int rightface[12]   = {L,  T,  N,  L,  B,  R,  R,  F,  B,  F,  N,  T};
  561.                  /* face on right when going corner1 to corner2 */
  562.  
  563.  
  564. /* docube: triangulate the cube directly, without decomposition */
  565.  
  566. int docube (cube, p)
  567. CUBE *cube;
  568. PROCESS *p;
  569. {
  570.     INTLISTS *polys;
  571.     int i, index = 0;
  572.     for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0) index += (1<<i);
  573.     for (polys = cubetable[index]; polys; polys = polys->next) {
  574.     INTLIST *edges;
  575.     int a = -1, b = -1, count = 0;
  576.     for (edges = polys->list; edges; edges = edges->next) {
  577.         CORNER *c1 = cube->corners[corner1[edges->i]];
  578.         CORNER *c2 = cube->corners[corner2[edges->i]];
  579.         int c = vertid(c1, c2, p);
  580.         if (++count > 2 && ! p->triproc(a, b, c, p->vertices)) return 0;
  581.         if (count < 3) a = b;
  582.         b = c;
  583.     }
  584.     }
  585.     return 1;
  586. }
  587.  
  588.  
  589. /* nextcwedge: return next clockwise edge from given edge around given face */
  590.  
  591. int nextcwedge (edge, face)
  592. int edge, face;
  593. {
  594.     switch (edge) {
  595.     case LB: return (face == L)? LF : BN;
  596.     case LT: return (face == L)? LN : TF;
  597.     case LN: return (face == L)? LB : TN;
  598.     case LF: return (face == L)? LT : BF;
  599.     case RB: return (face == R)? RN : BF;
  600.     case RT: return (face == R)? RF : TN;
  601.     case RN: return (face == R)? RT : BN;
  602.     case RF: return (face == R)? RB : TF;
  603.     case BN: return (face == B)? RB : LN;
  604.     case BF: return (face == B)? LB : RF;
  605.     case TN: return (face == T)? LT : RN;
  606.     case TF: return (face == T)? RT : LF;
  607.     }
  608. }
  609.  
  610.  
  611. /* otherface: return face adjoining edge that is not the given face */
  612.  
  613. int otherface (edge, face)
  614. int edge, face;
  615. {
  616.     int other = leftface[edge];
  617.     return face == other? rightface[edge] : other;
  618. }
  619.  
  620.  
  621. /* makecubetable: create the 256 entry table for cubical polygonization */
  622.  
  623. makecubetable ()
  624. {
  625.     int i, e, c, done[12], pos[8];
  626.     for (i = 0; i < 256; i++) {
  627.     for (e = 0; e < 12; e++) done[e] = 0;
  628.     for (c = 0; c < 8; c++) pos[c] = BIT(i, c);
  629.     for (e = 0; e < 12; e++)
  630.         if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
  631.         INTLIST *ints = 0;
  632.         INTLISTS *lists = (INTLISTS *) mycalloc(1, sizeof(INTLISTS));
  633.         int start = e, edge = e;
  634.         /* get face that is to right of edge from pos to neg corner: */
  635.         int face = pos[corner1[e]]? rightface[e] : leftface[e];
  636.         while (1) {
  637.             edge = nextcwedge(edge, face);
  638.             done[edge] = 1;
  639.             if (pos[corner1[edge]] != pos[corner2[edge]]) {
  640.             INTLIST *tmp = ints;
  641.             ints = (INTLIST *) mycalloc(1, sizeof(INTLIST));
  642.             ints->i = edge;
  643.             ints->next = tmp; /* add edge to head of list */
  644.             if (edge == start) break;
  645.             face = otherface(edge, face);
  646.             }
  647.         }
  648.         lists->list = ints; /* add ints to head of table entry */
  649.         lists->next = cubetable[i];
  650.         cubetable[i] = lists;
  651.         }
  652.     }
  653. }
  654.  
  655.  
  656. /**** Storage ****/
  657.  
  658.  
  659. /* mycalloc: return successful calloc or exit program */
  660.  
  661. char *mycalloc (nitems, nbytes)
  662. int nitems, nbytes;
  663. {
  664.    char *ptr = calloc(nitems, nbytes);
  665.    if (ptr != NULL) return ptr;
  666.    fprintf(stderr, "can't calloc %d bytes\n", nitems*nbytes);
  667.    exit(1);
  668. }
  669.  
  670.  
  671. /* setcenter: set (i,j,k) entry of table[]
  672.  * return 1 if already set; otherwise, set and return 0 */
  673.  
  674. int setcenter(table, i, j, k)
  675. CENTERLIST *table[];
  676. int i, j, k;
  677. {
  678.     int index = HASH(i, j, k);
  679.     CENTERLIST *new, *l, *q = table[index];
  680.     for (l = q; l != NULL; l = l->next)
  681.     if (l->i == i && l->j == j && l->k == k) return 1;
  682.     new = (CENTERLIST *) mycalloc(1, sizeof(CENTERLIST));
  683.     new->i = i; new->j = j; new->k = k; new->next = q;
  684.     table[index] = new;
  685.     return 0;
  686. }
  687.  
  688.  
  689. /* setedge: set vertex id for edge */
  690.  
  691. setedge (table, i1, j1, k1, i2, j2, k2, vid)
  692. EDGELIST *table[];
  693. int i1, j1, k1, i2, j2, k2, vid;
  694. {
  695.     unsigned int index;
  696.     EDGELIST *new;
  697.     if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
  698.     int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t;
  699.     }
  700.     index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
  701.     new = (EDGELIST *) mycalloc(1, sizeof(EDGELIST));
  702.     new->i1 = i1; new->j1 = j1; new->k1 = k1;
  703.     new->i2 = i2; new->j2 = j2; new->k2 = k2;
  704.     new->vid = vid;
  705.     new->next = table[index];
  706.     table[index] = new;
  707. }
  708.  
  709.  
  710. /* getedge: return vertex id for edge; return -1 if not set */
  711.  
  712. int getedge (table, i1, j1, k1, i2, j2, k2)
  713. EDGELIST *table[];
  714. int i1, j1, k1, i2, j2, k2;
  715. {
  716.     EDGELIST *q;
  717.     if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
  718.     int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t;
  719.     };
  720.     q = table[HASH(i1, j1, k1)+HASH(i2, j2, k2)];
  721.     for (; q != NULL; q = q->next)
  722.     if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
  723.         q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
  724.         return q->vid;
  725.     return -1;
  726. }
  727.  
  728.  
  729. /**** Vertices ****/
  730.  
  731.  
  732. /* vertid: return index for vertex on edge:
  733.  * c1->value and c2->value are presumed of different sign
  734.  * return saved index if any; else compute vertex and save */
  735.  
  736. int vertid (c1, c2, p)
  737. CORNER *c1, *c2;
  738. PROCESS *p;
  739. {
  740.     VERTEX v;
  741.     POINT a, b;
  742.     int vid = getedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
  743.     if (vid != -1) return vid;                 /* previously computed */
  744.     a.x = c1->x; a.y = c1->y; a.z = c1->z;
  745.     b.x = c2->x; b.y = c2->y; b.z = c2->z;
  746.     converge(&a, &b, c1->value, p->function, &v.position); /* position */
  747.     vnormal(&v.position, p, &v.normal);               /* normal */
  748.     addtovertices(&p->vertices, v);               /* save vertex */
  749.     vid = p->vertices.count-1;
  750.     setedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
  751.     return vid;
  752. }
  753.  
  754.  
  755. /* addtovertices: add v to sequence of vertices */
  756.  
  757. addtovertices (vertices, v)
  758. VERTICES *vertices;
  759. VERTEX v;
  760. {
  761.     if (vertices->count == vertices->max) {
  762.     int i;
  763.     VERTEX *new;
  764.     vertices->max = vertices->count == 0 ? 10 : 2*vertices->count;
  765.     new = (VERTEX *) mycalloc(vertices->max, sizeof(VERTEX));
  766.     for (i = 0; i < vertices->count; i++) new[i] = vertices->ptr[i];
  767.     if (vertices->ptr != NULL) free((char *) vertices->ptr);
  768.     vertices->ptr = new;
  769.     }
  770.     vertices->ptr[vertices->count++] = v;
  771. }
  772.  
  773.  
  774. /* vnormal: compute unit length surface normal at point */
  775.  
  776. vnormal (point, p, v)
  777. POINT *point, *v;
  778. PROCESS *p;
  779. {
  780.     double f = p->function(point->x, point->y, point->z);
  781.     v->x = p->function(point->x+p->delta, point->y, point->z)-f;
  782.     v->y = p->function(point->x, point->y+p->delta, point->z)-f;
  783.     v->z = p->function(point->x, point->y, point->z+p->delta)-f;
  784.     f = sqrt(v->x*v->x + v->y*v->y + v->z*v->z);
  785.     if (f != 0.0) {v->x /= f; v->y /= f; v->z /= f;}
  786. }
  787.  
  788.  
  789. /* converge: from two points of differing sign, converge to zero crossing */
  790.  
  791. converge (p1, p2, v, function, p)
  792. double v;
  793. double (*function)();
  794. POINT *p1, *p2, *p;
  795. {
  796.     int i = 0;
  797.     POINT pos, neg;
  798.     if (v < 0) {
  799.     pos.x = p2->x; pos.y = p2->y; pos.z = p2->z;
  800.     neg.x = p1->x; neg.y = p1->y; neg.z = p1->z;
  801.     }
  802.     else {
  803.     pos.x = p1->x; pos.y = p1->y; pos.z = p1->z;
  804.     neg.x = p2->x; neg.y = p2->y; neg.z = p2->z;
  805.     }
  806.     while (1) {
  807.     p->x = 0.5*(pos.x + neg.x);
  808.     p->y = 0.5*(pos.y + neg.y);
  809.     p->z = 0.5*(pos.z + neg.z);
  810.     if (i++ == RES) return;
  811.     if ((function(p->x, p->y, p->z)) > 0.0)
  812.          {pos.x = p->x; pos.y = p->y; pos.z = p->z;}
  813.     else {neg.x = p->x; neg.y = p->y; neg.z = p->z;}
  814.     }
  815. }
  816.